python3.11 -m venv hypoxia_env

# Loading Libraries
import os # for file path manipulations
import glob # for file pattern matching

# Third-party Libraries
import matplotlib.pyplot as plt # for plotting
import seaborn as sns # for statistical data visualization
from matplotlib.colors import LogNorm # for logarithmic color normalization
import pandas as pd # for data manipulation
import numpy as np # for numerical operations
from scipy import stats # for statistical functions

import xarray as xr  # for working with labeled multi-dimensional arrays
import rasterio # for raster data processing
import ee # for Google Earth Engine
import geemap # for geospatial data analysis and visualization

# Set visualization defaults
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (12, 6)
np.random.seed(42)

df = pd.read_csv(r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\NC\Dissolved_Oxygen_-_Gulf_of_Mexico_(GCOOS).csv", usecols=['resultmeasurevalue_mg_l','latitude','longitude','activitystartdate'],parse_dates=['activitystartdate'])
df.head(3)

# select half of rows if 'resultmeasurevalue_mg_l' <2 and half of rows if 'resultmeasurevalue_mg_l' >2
# remove rows if 'resultmeasurevalue_mg_l' > 10
df2=df[df['resultmeasurevalue_mg_l']<10]
# greater than 0 : negative DO is not needed
df2=df2[df2['resultmeasurevalue_mg_l']>0]
# remove rows if -99, -9999, -999, -9.999, -999.999
df2=df2[~df2['resultmeasurevalue_mg_l'].isin([-99,-9999,-999,-9.999,-999.999])]

# Exploratory Data Analysis: Distribution and Statistics
# Basic statistics
print("=" * 60)
print("BASIC STATISTICS")
print("=" * 60)
print(df.describe())
print("\n" + "=" * 60)
print("DATA TYPES & MISSING VALUES")
print("=" * 60)
print(df.info())
print("\n" + "=" * 60)
print("MISSING VALUES COUNT")
print("=" * 60)
print(df.isnull().sum())

============================================================
BASIC STATISTICS
============================================================
           latitude     longitude  resultmeasurevalue_mg_l          year
count  6.725953e+06  6.725953e+06             6.725953e+06  6.725953e+06
mean   2.885543e+01 -8.780485e+01            -1.479944e+03  2.007931e+03
std    1.721641e+00  4.956432e+00             1.219334e+06  5.130250e+00
min    1.879200e+01 -9.772330e+01            -1.000000e+09  1.922000e+03
25%    2.794560e+01 -8.846290e+01             4.900000e+00  2.007000e+03
50%    2.970210e+01 -8.782280e+01             6.500000e+00  2.009000e+03
75%    3.038080e+01 -8.487520e+01             8.100000e+00  2.011000e+03
max    3.084060e+01 -8.000000e+01             9.999900e+04  2.014000e+03

============================================================
DATA TYPES & MISSING VALUES
============================================================
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6725953 entries, 0 to 6725952
Data columns (total 5 columns):
 #   Column                   Dtype              
---  ------                   -----              
 0   latitude                 float64            
 1   longitude                float64            
 2   activitystartdate        datetime64[ns, UTC]
 3   resultmeasurevalue_mg_l  float64            
 4   year                     int32              
dtypes: datetime64[ns, UTC](1), float64(3), int32(1)
memory usage: 230.9 MB
None

============================================================
MISSING VALUES COUNT
============================================================
latitude                   0
longitude                  0
activitystartdate          0
resultmeasurevalue_mg_l    0
year                       0
dtype: int64

# Distribution of dissolved oxygen values
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Histogram of DO values
axes[0, 0].hist(df2['resultmeasurevalue_mg_l'].dropna(), bins=50, color='steelblue', edgecolor='black')
axes[0, 0].axvline(2, color='red', linestyle='--', linewidth=2, label='Hypoxia Threshold (2 mg/L)')
axes[0, 0].set_xlabel('Dissolved Oxygen (mg/L)', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('Frequency', fontsize=12, fontweight='bold')
axes[0, 0].set_title('Distribution of Dissolved Oxygen Values', fontsize=14, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)

# Box plot
axes[0, 1].boxplot(df2['resultmeasurevalue_mg_l'].dropna(), vert=True)
axes[0, 1].axhline(2, color='red', linestyle='--', linewidth=2, label='Hypoxia Threshold')
axes[0, 1].set_ylabel('Dissolved Oxygen (mg/L)', fontsize=12, fontweight='bold')
axes[0, 1].set_title('DO Value Box Plot (Outliers Visible)', fontsize=14, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)

# Spatial distribution scatter
scatter = axes[1, 0].scatter(df['longitude'], df['latitude'], 
                             c=df['resultmeasurevalue_mg_l'], 
                             cmap='RdYlGn', s=1, alpha=0.6, vmin=0, vmax=10)
axes[1, 0].set_xlabel('Longitude', fontsize=12, fontweight='bold')
axes[1, 0].set_ylabel('Latitude', fontsize=12, fontweight='bold')
axes[1, 0].set_title('Spatial Distribution of DO Measurements', fontsize=14, fontweight='bold')
plt.colorbar(scatter, ax=axes[1, 0], label='DO (mg/L)')
axes[1, 0].grid(alpha=0.3)

# Temporal distribution
df['year'] = df['activitystartdate'].dt.year
year_counts = df['year'].value_counts().sort_index()
axes[1, 1].plot(year_counts.index, year_counts.values, marker='o', linewidth=2, markersize=5, color='darkgreen')
axes[1, 1].set_xlabel('Year', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('Number of Measurements', fontsize=12, fontweight='bold')
axes[1, 1].set_title('Temporal Coverage of DO Measurements', fontsize=14, fontweight='bold')
axes[1, 1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# Hypoxia statistics
hypoxic_count = (df['resultmeasurevalue_mg_l'] < 2).sum()
total_count = df['resultmeasurevalue_mg_l'].notna().sum()
hypoxic_pct = (hypoxic_count / total_count * 100) if total_count > 0 else 0

print("\n" + "=" * 60)
print("HYPOXIA STATISTICS")
print("=" * 60)
print(f"Total measurements: {total_count:,}")
print(f"Hypoxic measurements (DO < 2 mg/L): {hypoxic_count:,}")
print(f"Hypoxic percentage: {hypoxic_pct:.2f}%")
print(f"Date range: {df['activitystartdate'].min()} to {df['activitystartdate'].max()}")
print(f"Spatial extent: Lat [{df['latitude'].min():.2f}, {df['latitude'].max():.2f}], Lon [{df['longitude'].min():.2f}, {df['longitude'].max():.2f}]")
print("=" * 60)

============================================================
HYPOXIA STATISTICS
============================================================
Total measurements: 6,725,953
Hypoxic measurements (DO < 2 mg/L): 296,657
Hypoxic percentage: 4.41%
Date range: 1922-02-01 00:00:00+00:00 to 2014-08-14 00:00:00+00:00
Spatial extent: Lat [18.79, 30.84], Lon [-97.72, -80.00]
============================================================

# Add binary column for hypoxia classification
# Creates a binary label: 1 = non-hypoxic (DO > 2 mg/L), 0 = hypoxic (DO ≤ 2 mg/L)
# Note: This encoding is later inverted during model training where we use (DO < 2) as class 1 (hypoxic)
df2['hypoxia'] = (df2['resultmeasurevalue_mg_l'] > 2).astype(int)

# Convert 'activitystartdate' to datetime
df2['activitystartdate'] = pd.to_datetime(df2['activitystartdate'])

# Extract year from 'activitystartdate'
df2['year'] = df2['activitystartdate'].dt.year
df2.head(3)

df2=df2[df2['year']>1970]
# Count hypoxia events per year #To count per month: df2.groupby(df2['activitystartdate'].dt.to_period("M"))['hypoxia'].sum().reset_index()
hypoxia_counts = df2.groupby('year')['hypoxia'].sum().reset_index()

# Plotting
plt.figure(figsize=(12, 6))
sns.lineplot(x='year', y='hypoxia', data=hypoxia_counts, marker='o', label='Hypoxia Events')

# Adding regression line
sns.regplot(x='year', y='hypoxia', data=hypoxia_counts, scatter=False, label='Regression Line')
# Set y-axis to log scale
plt.title('Hypoxia Events Per Year with Regression Line', fontsize=14,fontweight='bold')
plt.xlabel('Year', fontsize=12,fontweight='bold')
plt.ylabel('Number of Hypoxia Events Per Years (Combined)', fontsize=12,fontweight='bold')
# xlim minimum and maximum of df2['year']
plt.xlim(df2['year'].min()-1, df2['year'].max()+2)
# ylimit 0 to maximum of hypoxia_counts['hypoxia'] +10
plt.ylim(0, hypoxia_counts['hypoxia'].max()+10)
# add xticks
plt.xticks(np.arange(df2['year'].min()-1, df2['year'].max()+2, 3.0))
plt.legend()
# save figure 
#plt.savefig(r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\figure\hypoxiaevtn.jpg",dpi=400,bbox_inches='tight')
plt.show()

# Create a new column for the 5/10-year periods
df2['5_year_period'] = (df2['activitystartdate'].dt.year // 10) * 10

# Group by the 5-year periods and calculate descriptive statistics
hypoxia_stats = df2.groupby('5_year_period')['hypoxia'].describe()

print(hypoxia_stats)

                   count      mean       std  min  25%  50%  75%  max
5_year_period                                                        
1970              8102.0  0.997902  0.045761  0.0  1.0  1.0  1.0  1.0
1980             26517.0  0.997322  0.051676  0.0  1.0  1.0  1.0  1.0
1990            281893.0  0.977438  0.148502  0.0  1.0  1.0  1.0  1.0
2000           3119826.0  0.954849  0.207636  0.0  1.0  1.0  1.0  1.0
2010           2656997.0  0.949562  0.218847  0.0  1.0  1.0  1.0  1.0

study_area = ee.FeatureCollection("projects/ee-hafezahmad10/assets/US_GULF_eez_FL_MS_L_T")
import calendar
years = range(2018, 2022)  # Define the range of years
date_start = []
date_end = []
for year in years:
    for month in range(1, 13):  # Loop over each month (1 to 12)
        # Get the number of days in the month
        num_days = calendar.monthrange(year, month)[1]
        # Append the start and end dates for the current month
        date_start.append(f'{year}-{month:02d}-01')
        date_end.append(f'{year}-{month:02d}-{num_days}')
# now use for loop to iterate over date_start and date_end and image_collection
for start, end in zip(date_start, date_end):
    # filter image collection
    image_collection = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI').filterBounds(study_area).filterDate(start, end).mean().clip(study_area)
    # generate text B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 and B8A
    bands=['chlor_a', 'nflh','poc','sst','Rrs_412','Rrs_443','Rrs_469','Rrs_488','Rrs_531','Rrs_547','Rrs_555','Rrs_645','Rrs_667','Rrs_678']
    out_tif =f"OCEANDATA_MODIS-Aqua_L3SMI{start}"
    # to drive
    geemap.ee_export_image_to_drive(image_collection.select(bands),description=out_tif, scale=1000, region=study_area.geometry(),folder='GEE_exports')

# add new column:id:sequence of numbers starting from 1
df2['id']=np.arange(1,len(df2)+1)
df2

# Convert DataFrame to Earth Engine FeatureCollection
point2022=geemap.df_to_ee(df2)
study_area = ee.FeatureCollection("projects/ee-hafezahmad10/assets/US_GULF_eez_FL_MS_L_T")
image_collection = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI').filterBounds(study_area).filterDate('2010-01-01', '2010-01-31').mean().clip(study_area)
# generate text B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 and B8A
bands=['chlor_a', 'nflh','poc','sst','Rrs_412','Rrs_443','Rrs_469','Rrs_488','Rrs_531','Rrs_547','Rrs_555','Rrs_645','Rrs_667','Rrs_678']
path=r"D:\my works\bloom\data\MODIS"
out_csv = os.path.join(path, 'train.csv')
# downloading data to csv ,scale 500 meters
geemap.extract_values_to_points(point2022, image_collection.select(bands), out_csv,scale=500)

path=r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\csv"
# find all "train" containing .csv files
files = glob.glob(os.path.join(path, "*train*.csv"))
# read first only
df=pd.read_csv(files[0])
df.drop(columns=['OID_', 'latitude', 'longitude', 'activitystartdate'], inplace=True)
# read next all and rename columns based on first file columns
for file in files[1:]:
    df1 = pd.read_csv(file)
    df1.drop(columns=['OID_', 'latitude', 'longitude', 'activitystartdate'], inplace=True)
    df1.columns = df.columns
    df = pd.concat([df, df1], ignore_index=True)
# remove nan values
df.dropna(inplace=True)
print(df.shape)
print(df.columns)
df.head(2)

(116444, 15)
Index(['resultmeasurevalue_mg_l', 'chlor_a', 'nflh', 'poc', 'sst', 'Rrs_412',
       'Rrs_443', 'Rrs_469', 'Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555',
       'Rrs_645', 'Rrs_667', 'Rrs_678'],
      dtype='object')

# Model training and evaluation
import joblib # for saving models
from sklearn.model_selection import train_test_split # for splitting data
from sklearn.ensemble import RandomForestClassifier # for ensemble learning
from sklearn.linear_model import LogisticRegression # for logistic regression
from sklearn.svm import SVC # for support vector machine
from sklearn.neighbors import KNeighborsClassifier # 
from sklearn.tree import DecisionTreeClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.neural_network import MLPClassifier
from sklearn.metrics import (
    accuracy_score,
    balanced_accuracy_score,
    classification_report,
    precision_recall_fscore_support,
    roc_auc_score,
)
from sklearn.preprocessing import MinMaxScaler
from sklearn.pipeline import make_pipeline
from sklearn.utils import resample
from xgboost import XGBClassifier

# Reusable scaler (used only inside pipelines)
scaler = MinMaxScaler() # scaling features to a given range

# Feature and target definitions
feature_columns = [
    'chlor_a', 'nflh', 'poc', 'sst', 'Rrs_412', 'Rrs_443', 'Rrs_469',
    'Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555', 'Rrs_645', 'Rrs_667', 'Rrs_678'
]
target_column = 'resultmeasurevalue_mg_l'

# Binary label: 1 = hypoxic (<2 mg/L), 0 = non-hypoxic

def convert_to_class(value: float) -> int:
    return 1 if value < 2 else 0

df['class_label'] = df[target_column].apply(convert_to_class)
df.drop(columns=[target_column], inplace=True)
print(df['class_label'].value_counts())

class_label
0    115248
1      1196
Name: count, dtype: int64

# Address class imbalance by downsampling the majority class
majority = df[df.class_label == 0]# majority class: non-hypoxic
minority = df[df.class_label == 1] # minority class: hypoxic
majority_downsampled = resample(
    majority,
    replace=False,
    n_samples=minority.shape[0],
    random_state=123,
)
df_balanced = pd.concat([minority, majority_downsampled])
print(df_balanced.class_label.value_counts())

class_label
1    1196
0    1196
Name: count, dtype: int64

# Split features/labels
X = df_balanced[feature_columns]
y = df_balanced['class_label']

# Train/test split with stratification
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
) #  test size 30% and train size 70%

# Create and save a standalone MinMaxScaler fitted on training data
# This scaler will be used for feature normalization during prediction on new raster data
from sklearn.preprocessing import MinMaxScaler

scaler = MinMaxScaler()
scaler.fit(X_train)
models_directory = r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1'
# Persist the fitted scaler
scaler_path = os.path.join(models_directory, 'MinMaxScaler.pkl')
joblib.dump(scaler, scaler_path)

print(f'✓ Created and saved MinMaxScaler to: {scaler_path}')
print(f'  Scaler fitted on training data with shape: {X_train.shape}')
print(f'  Feature min/max ranges: {list(zip(scaler.data_min_[:3], scaler.data_max_[:3]))}...')

# Model registry with simple, sane defaults
models = {
    'Random Forest': {
        'est': RandomForestClassifier(
            n_estimators=300, # n_estimators number of trees
            max_depth=None,
            n_jobs=-1, # use all processors
            class_weight='balanced',
            random_state=42,# really important for reproducibility
        ),
        'scale': False,
    },
    'Logistic Regression': {
        'est': LogisticRegression(
            max_iter=1000, # increase max iterations
            class_weight='balanced',
            solver='lbfgs',# good for small datasets,lbfgs stands for Limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm
            random_state=42,
        ),
        'scale': True,
    },
    'Support Vector Machine': {
        'est': SVC(
            kernel='rbf',# radial basis function kernel, which is good for non-linear data
            C=2.0,
            gamma='scale',# scale is 1 / (n_features * X.var()) as value of gamma
            probability=True,
            class_weight='balanced',# important for imbalanced datasets
            random_state=42,
        ),
        'scale': True,
    },
    'K-Nearest Neighbors': {
        'est': KNeighborsClassifier(n_neighbors=7, weights='distance'),# weights='distance' means closer neighbors have more influence
        'scale': True,
    },
    'Decision Tree': {
        'est': DecisionTreeClassifier(
            max_depth=8,# limit depth to prevent overfitting
            class_weight='balanced',
            random_state=42,
        ),
        'scale': False,
    },
    'Naive Bayes': {
        'est': GaussianNB(),# Gaussian Naive Bayes classifier
        'scale': True,
    },
    'XGBoost': {
        'est': XGBClassifier(
            n_estimators=300,
            max_depth=4,
            learning_rate=0.05,# smaller learning rate for better convergence
            subsample=0.8, # subsample ratio of the training instances
            colsample_bytree=0.8,
            eval_metric='logloss',# evaluation metric
            seed=42,
        ),
        'scale': False,
    },
    'MLP': {
        'est': MLPClassifier(
            hidden_layer_sizes=(64, 32),# 64, 32 means we are using two hidden layers with 64 and 32 neurons respectively
            max_iter=500, # maximum number of iterations
            random_state=42,
        ),
        'scale': True,
    },
}
# models directory
models_directory = r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1'
os.makedirs(models_directory, exist_ok=True)

# Function to compute evaluation metrics: accuracy, balanced accuracy, precision, recall, F1 score, and ROC AUC
def get_scores(model, X_te, y_te):
    y_pred = model.predict(X_te)
    proba = None
    if hasattr(model, 'predict_proba'):
        proba = model.predict_proba(X_te)[:, 1]
    elif hasattr(model, 'decision_function'):
        proba = model.decision_function(X_te)

    acc = accuracy_score(y_te, y_pred)
    bal_acc = balanced_accuracy_score(y_te, y_pred)
    precision, recall, f1, _ = precision_recall_fscore_support(
        y_te, y_pred, average='binary', zero_division=0
    )
    roc = roc_auc_score(y_te, proba) if proba is not None else np.nan
    return acc, bal_acc, precision, recall, f1, roc

results = []
# Train, evaluate, and persist each model
for name, cfg in models.items():
    est = cfg['est']
    model = make_pipeline(MinMaxScaler(), est) if cfg['scale'] else est
    model.fit(X_train, y_train)

    acc, bal_acc, precision, recall, f1, roc = get_scores(model, X_test, y_test)

    # Persist model
    model_path = os.path.join(models_directory, f'{name}.pkl')
    joblib.dump(model, model_path)

    results.append({
        'Classifier': name,
        'Accuracy': acc,
        'Balanced_Accuracy': bal_acc,
        'Precision': precision,
        'Recall': recall,
        'F1': f1,
        'ROC_AUC': roc,
        'Model_Path': model_path,
    })
# Create results DataFrame and sort by F1 score
results_df = pd.DataFrame(results)
# show results sorted by F1 score
results_df.sort_values(by='F1', ascending=False)

# Optional: hyperparameter tuning for Random Forest
from sklearn.model_selection import RandomizedSearchCV

rf_base = RandomForestClassifier(
    n_estimators=400,
    class_weight='balanced',
    random_state=42,
    n_jobs=-1,
)
param_dist = {
    'n_estimators': [200, 300, 400, 500],
    'max_depth': [None, 6, 8, 10, 12],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4],
    'max_features': ['sqrt', 'log2'],
}
# Perform randomized search
search = RandomizedSearchCV(
    rf_base,
    param_distributions=param_dist,
    n_iter=12,
    cv=3,
    scoring='f1',  # prioritize hypoxia class detection
    n_jobs=-1,
    random_state=42,
    verbose=0,
)
# Perform randomized search
search.fit(X_train, y_train)

best_rf = search.best_estimator_
acc, bal_acc, precision, recall, f1, roc = get_scores(best_rf, X_test, y_test)

# Persist tuned model
rf_tuned_path = os.path.join(models_directory, 'RandomForest_tuned.pkl')
joblib.dump(best_rf, rf_tuned_path)

# Append tuned results to the summary table
rf_row = {
    'Classifier': 'Random Forest (tuned)',
    'Accuracy': acc,
    'Balanced_Accuracy': bal_acc,
    'Precision': precision,
    'Recall': recall,
    'F1': f1,
    'ROC_AUC': roc,
    'Model_Path': rf_tuned_path,
}
results_df = pd.concat([results_df, pd.DataFrame([rf_row])], ignore_index=True)

# Choose overall best by F1 and save a canonical best_model.pkl
best_row_overall = results_df.sort_values(by='F1', ascending=False).iloc[0]
best_model_path_overall = best_row_overall['Model_Path']
best_model_overall = joblib.load(best_model_path_overall)
best_model_canonical_path = os.path.join(models_directory, 'best_model.pkl')
joblib.dump(best_model_overall, best_model_canonical_path)

print('Best RF params:', search.best_params_)
print('Tuned RF test metrics:', rf_row)
print('Overall best model:', best_row_overall['Classifier'])
print('Saved canonical best model to:', best_model_canonical_path)

results_df.sort_values(by='F1', ascending=False)

Best RF params: {'n_estimators': 200, 'min_samples_split': 10, 'min_samples_leaf': 1, 'max_features': 'sqrt', 'max_depth': 10}
Tuned RF test metrics: {'Classifier': 'Random Forest (tuned)', 'Accuracy': 0.915041782729805, 'Balanced_Accuracy': 0.915041782729805, 'Precision': 0.8860103626943006, 'Recall': 0.9526462395543176, 'F1': 0.9181208053691275, 'ROC_AUC': 0.9583996089415817, 'Model_Path': 'C:\\Users\\hafez\\MSU\\Research\\msGOM\\mssound\\bloom\\data\\models1\\RandomForest_tuned.pkl'}
Overall best model: Random Forest
Saved canonical best model to: C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1\best_model.pkl

# Visualize model metrics and inspect best model
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import ConfusionMatrixDisplay, RocCurveDisplay, classification_report

# Bar plot of metrics across models
metrics_to_plot = ['Accuracy', 'Balanced_Accuracy', 'Precision', 'Recall', 'F1', 'ROC_AUC']
plot_df = results_df.melt(id_vars=['Classifier'], value_vars=metrics_to_plot,
                          var_name='Metric', value_name='Score')
plt.figure(figsize=(12, 6))
sns.barplot(data=plot_df, x='Classifier', y='Score', hue='Metric')
plt.xticks(rotation=30, ha='right')
plt.ylim(0, 1)
plt.title('Model Performance Comparison')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

# Choose best model by F1
best_row = results_df.sort_values(by='F1', ascending=False).iloc[0]
best_model_path = best_row['Model_Path']
best_model_name = best_row['Classifier']
best_model = joblib.load(best_model_path)

# Predictions for diagnostics
y_pred = best_model.predict(X_test)

# Classification report
print(f"Best model: {best_model_name}\n")
print(classification_report(y_test, y_pred, digits=3, zero_division=0))

# Confusion matrix
ConfusionMatrixDisplay.from_predictions(y_test, y_pred, display_labels=['Non-hypoxic','Hypoxic'], cmap='Blues')
plt.title(f'Confusion Matrix: {best_model_name}')
plt.tight_layout()
plt.show()

# ROC curve (skip if no scores available)
y_scores = None
if hasattr(best_model, 'predict_proba'):
    y_scores = best_model.predict_proba(X_test)[:, 1]
elif hasattr(best_model, 'decision_function'):
    y_scores = best_model.decision_function(X_test)

if y_scores is not None:
    RocCurveDisplay.from_predictions(y_test, y_scores, name=best_model_name)
    plt.title(f'ROC Curve: {best_model_name}')
    plt.tight_layout()
    plt.show()
else:
    print('ROC curve skipped: model does not provide scores/probabilities.')

Best model: Random Forest

              precision    recall  f1-score   support

           0      0.949     0.877     0.912       359
           1      0.886     0.953     0.918       359

    accuracy                          0.915       718
   macro avg      0.917     0.915     0.915       718
weighted avg      0.917     0.915     0.915       718

# Bootstrap uncertainty on test metrics for the best model
from sklearn.metrics import precision_score, recall_score, f1_score, roc_auc_score

# Requires best_model, X_test, y_test from previous cell
n_boot = 200
rng = np.random.default_rng(42)

has_scores = hasattr(best_model, 'predict_proba') or hasattr(best_model, 'decision_function')

boot_metrics = []
for _ in range(n_boot):
    idx = rng.choice(len(y_test), size=len(y_test), replace=True)
    Xb = X_test.iloc[idx] if hasattr(X_test, 'iloc') else X_test[idx]
    yb = y_test.iloc[idx] if hasattr(y_test, 'iloc') else y_test[idx]

    y_pred_b = best_model.predict(Xb)
    if hasattr(best_model, 'predict_proba'):
        y_score_b = best_model.predict_proba(Xb)[:, 1]
    elif hasattr(best_model, 'decision_function'):
        y_score_b = best_model.decision_function(Xb)
    else:
        y_score_b = None

    prec = precision_score(yb, y_pred_b, zero_division=0)
    rec = recall_score(yb, y_pred_b, zero_division=0)
    f1 = f1_score(yb, y_pred_b, zero_division=0)
    roc = roc_auc_score(yb, y_score_b) if y_score_b is not None else np.nan
    boot_metrics.append([prec, rec, f1, roc])

boot_arr = np.array(boot_metrics)
quantiles = np.nanpercentile(boot_arr, [2.5, 50, 97.5], axis=0)

metric_names = ['Precision', 'Recall', 'F1', 'ROC_AUC']
print('Bootstrap 95% CIs (test set):')
for i, name in enumerate(metric_names):
    lo, med, hi = quantiles[:, i]
    print(f"{name}: median={med:.3f}, 95% CI=({lo:.3f}, {hi:.3f})")

Bootstrap 95% CIs (test set):
Precision: median=0.885, 95% CI=(0.853, 0.916)
Recall: median=0.953, 95% CI=(0.931, 0.971)
F1: median=0.918, 95% CI=(0.894, 0.936)
ROC_AUC: median=0.962, 95% CI=(0.946, 0.973)

# Feature importance and SHAP for the best model


# Helper to grab the final estimator from a pipeline or standalone model
def get_final_estimator(model):
    if hasattr(model, 'named_steps'):
        # walk reversed to get last estimator-like step
        for step in reversed(list(model.named_steps.values())):
            if hasattr(step, 'feature_importances_') or hasattr(step, 'coef_'):
                return step
        return list(model.named_steps.values())[-1]
    return model

final_est = get_final_estimator(best_model)

# Build a DataFrame of feature values (for plotting and SHAP)
if hasattr(X_test, 'columns'):
    X_test_df = X_test.copy()
else:
    X_test_df = pd.DataFrame(X_test, columns=feature_columns)

# Permutation/embedded feature importance
importance_df = None
if hasattr(final_est, 'feature_importances_'):
    vals = final_est.feature_importances_
    importance_df = pd.DataFrame({'feature': feature_columns, 'importance': vals})
elif hasattr(final_est, 'coef_'):
    vals = np.abs(np.ravel(final_est.coef_))
    importance_df = pd.DataFrame({'feature': feature_columns, 'importance': vals})

if importance_df is not None:
    importance_df = importance_df.sort_values('importance', ascending=True)
    plt.figure(figsize=(8, 6))
    plt.barh(importance_df['feature'], importance_df['importance'])
    plt.title(f'Feature importance: {best_model_name}')
    plt.tight_layout()
    plt.show()
else:
    print('Feature importance not available for this model type.')

# Precision-Recall curve, PR-AUC, and threshold sweep for the best model
from sklearn.metrics import precision_recall_curve, average_precision_score

# Ensure we have probability scores
if hasattr(best_model, 'predict_proba'):
    y_score = best_model.predict_proba(X_test)[:, 1]
elif hasattr(best_model, 'decision_function'):
    y_score = best_model.decision_function(X_test)
else:
    raise ValueError('Best model does not expose decision scores or probabilities.')

# PR curve and PR-AUC (Average Precision)
precision_vals, recall_vals, thresh = precision_recall_curve(y_test, y_score)
pr_auc = average_precision_score(y_test, y_score)

# Sweep thresholds to maximize F1
f1_vals = 2 * (precision_vals * recall_vals) / (precision_vals + recall_vals + 1e-9)
best_idx = f1_vals.argmax()
best_threshold = thresh[best_idx - 1] if best_idx > 0 else thresh[0]

# Metrics at the chosen threshold
from sklearn.metrics import precision_score, recall_score, f1_score

y_pred_thresh = (y_score >= best_threshold).astype(int)
prec_t = precision_score(y_test, y_pred_thresh, zero_division=0)
rec_t = recall_score(y_test, y_pred_thresh, zero_division=0)
f1_t = f1_score(y_test, y_pred_thresh, zero_division=0)

print(f'PR-AUC (average precision): {pr_auc:.3f}')
print(f'Best threshold for F1: {best_threshold:.3f}')
print(f'Precision@thr: {prec_t:.3f}, Recall@thr: {rec_t:.3f}, F1@thr: {f1_t:.3f}')

# Plot PR curve with the chosen operating point
plt.figure(figsize=(8, 6))
plt.plot(recall_vals, precision_vals, label=f'PR curve (AP={pr_auc:.3f})')
plt.scatter(recall_vals[best_idx], precision_vals[best_idx], color='red', label='Chosen threshold')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.title(f'Precision-Recall: {best_model_name}')
plt.legend()
plt.tight_layout()
plt.show()

PR-AUC (average precision): 0.950
Best threshold for F1: 0.513
Precision@thr: 0.886, Recall@thr: 0.953, F1@thr: 0.918

# Reliability curve (calibration) and Brier score
from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.metrics import brier_score_loss
from sklearn.base import clone

# Base probabilities
if hasattr(best_model, 'predict_proba'):
    base_proba = best_model.predict_proba(X_test)[:, 1]
elif hasattr(best_model, 'decision_function'):
    base_proba = best_model.decision_function(X_test)
else:
    raise ValueError('Best model lacks probability/decision_function for calibration.')

brier_base = brier_score_loss(y_test, base_proba)

# Fit calibrated model (sigmoid/Platt scaling)
base_est = clone(best_model)
calibrator = CalibratedClassifierCV(base_est, method='sigmoid', cv=3)
calibrator.fit(X_train, y_train)
cal_proba = calibrator.predict_proba(X_test)[:, 1]
brier_cal = brier_score_loss(y_test, cal_proba)

print(f'Brier (uncalibrated): {brier_base:.4f}')
print(f'Brier (calibrated):   {brier_cal:.4f}')

# Calibration curve
prob_true_base, prob_pred_base = calibration_curve(y_test, base_proba, n_bins=10)
prob_true_cal, prob_pred_cal = calibration_curve(y_test, cal_proba, n_bins=10)

plt.figure(figsize=(7, 7))
plt.plot([0, 1], [0, 1], 'k--', label='Perfectly calibrated')
plt.plot(prob_pred_base, prob_true_base, marker='o', label='Uncalibrated')
plt.plot(prob_pred_cal, prob_true_cal, marker='o', label='Calibrated (sigmoid)')
plt.xlabel('Predicted probability')
plt.ylabel('Observed frequency')
plt.title(f'Reliability diagram: {best_model_name}')
plt.legend()
plt.tight_layout()
plt.show()

# Save calibrated model for downstream use
calibrated_model_path = os.path.join(models_directory, f'{best_model_name}_calibrated.pkl')
joblib.dump(calibrator, calibrated_model_path)
print('Saved calibrated model to:', calibrated_model_path)

Brier (uncalibrated): 0.0651
Brier (calibrated):   0.0681

Saved calibrated model to: C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1\Random Forest_calibrated.pkl

# Model-agnostic permutation importance and PDP/ICE plots
from sklearn.inspection import permutation_importance, PartialDependenceDisplay

# Permutation importance on the held-out test set
perm = permutation_importance(
    best_model,
    X_test,
    y_test,
    n_repeats=25,
    random_state=42,
    n_jobs=-1,
    scoring='f1',
)
perm_df = pd.DataFrame({
    'feature': feature_columns,
    'importance_mean': perm.importances_mean,
    'importance_std': perm.importances_std,
}).sort_values('importance_mean', ascending=True)

plt.figure(figsize=(8, 6))
plt.barh(perm_df['feature'], perm_df['importance_mean'], xerr=perm_df['importance_std'])
plt.title(f'Permutation importance (F1 loss): {best_model_name}')
plt.xlabel('Mean decrease in F1 when shuffled')
plt.tight_layout()
plt.show()

# Pick top features for PDP/ICE
top_features = perm_df.sort_values('importance_mean', ascending=False)['feature'].head(3).tolist()
print('Top features for PDP/ICE:', top_features)

fig, ax = plt.subplots(len(top_features), 1, figsize=(8, 4 * len(top_features)))
PartialDependenceDisplay.from_estimator(
    best_model,
    X_test,
    features=top_features,
    kind='both',  # PDP + ICE
    grid_resolution=30,
    ax=ax,
)
plt.tight_layout()
plt.show()

Top features for PDP/ICE: ['chlor_a', 'Rrs_488', 'Rrs_469']

# Partial Dependence and ICE for chlor_a
from sklearn.inspection import PartialDependenceDisplay

# Build a DataFrame with feature names if needed
if hasattr(X_test, 'columns'):
    X_pd = X_test
else:
    X_pd = pd.DataFrame(X_test, columns=feature_columns)

feature_name = 'chlor_a'
fig, ax = plt.subplots(1, 2, figsize=(12, 5))

# Average Partial Dependence
PartialDependenceDisplay.from_estimator(
    best_model, X_pd, [feature_name], ax=ax[0], kind='average', grid_resolution=50
)
ax[0].set_title('Partial Dependence: chlor_a')
ax[0].set_xlabel('chlor_a')
ax[0].set_ylabel('Predicted hypoxia probability')

# ICE (individual curves)
PartialDependenceDisplay.from_estimator(
    best_model, X_pd, [feature_name], ax=ax[1], kind='individual', subsample=200, grid_resolution=50
)
ax[1].set_title('ICE: chlor_a (subsample=200)')
ax[1].set_xlabel('chlor_a')
ax[1].set_ylabel('Predicted hypoxia probability')

plt.tight_layout()
plt.show()

# ROC curve for the best model
import os
import joblib
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, roc_auc_score

# Ensure X_test and y_test exist; ensure we have a model
try:
    best_model
except NameError:
    # Load canonical best model from disk
    models_directory = r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1'
    best_model_path = os.path.join(models_directory, 'best_model.pkl')
    best_model = joblib.load(best_model_path)

# Get scores for ROC (probabilities or decision function)
if hasattr(best_model, 'predict_proba'):
    y_score = best_model.predict_proba(X_test)[:, 1]
elif hasattr(best_model, 'decision_function'):
    y_score = best_model.decision_function(X_test)
    # Normalize decision function to [0, 1] range for plotting
    y_score = (y_score - y_score.min()) / (y_score.max() - y_score.min() + 1e-9)
else:
    raise ValueError('Best model does not expose probabilities or decision scores.')

fpr, tpr, _ = roc_curve(y_test, y_score)
auc_val = roc_auc_score(y_test, y_score)

plt.figure(figsize=(7, 6))
plt.plot(fpr, tpr, label=f'Best model ROC (AUC={auc_val:.3f})', color='#2563eb', linewidth=2)
plt.plot([0, 1], [0, 1], '--', color='gray', label='Chance')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate (1 - Specificity)')
plt.ylabel('True Positive Rate (Recall)')
plt.title('ROC Curve — Hypoxia Classification')
plt.legend(loc='lower right')
plt.grid(alpha=0.2)
plt.tight_layout()
plt.show()

# Overlay ROC curves for base RF, tuned RF, and calibrated model (if available)
import os
import joblib
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, roc_auc_score

# Require test data
if 'X_test' not in globals() or 'y_test' not in globals():
    print('⚠️ X_test/y_test not found. Please run the Train/Test Split cell first.')
else:
    models_directory = r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\models1'

    candidates = {
        'Base RF': os.path.join(models_directory, 'Random Forest.pkl'),
        'Tuned RF': os.path.join(models_directory, 'RandomForest_tuned.pkl'),
        'Calibrated': None,
    }
    # Try several calibrated filenames
    for cand in ['best_model_calibrated.pkl', 'RandomForest_tuned_calibrated.pkl', 'Random Forest_calibrated.pkl']:
        path_c = os.path.join(models_directory, cand)
        if os.path.exists(path_c):
            candidates['Calibrated'] = path_c
            break

    curves = {}
    for name, path in candidates.items():
        if path is None:
            print(f'ℹ️ {name} not found; skipping.')
            continue
        if not os.path.exists(path):
            print(f'ℹ️ {name} file missing: {path}; skipping.')
            continue
        try:
            model = joblib.load(path)
            if hasattr(model, 'predict_proba'):
                y_score = model.predict_proba(X_test)[:, 1]
            elif hasattr(model, 'decision_function'):
                y_score = model.decision_function(X_test)
                y_score = (y_score - y_score.min()) / (y_score.max() - y_score.min() + 1e-9)
            else:
                print(f'ℹ️ {name} has no probability/score method; skipping.')
                continue
            fpr, tpr, _ = roc_curve(y_test, y_score)
            auc_val = roc_auc_score(y_test, y_score)
            curves[name] = (fpr, tpr, auc_val)
        except Exception as e:
            print(f'⚠️ Could not evaluate {name}: {e}')

    # Plot
    if curves:
        plt.figure(figsize=(7.5, 6.5))
        palette = {
            'Base RF': '#374151',
            'Tuned RF': '#2563eb',
            'Calibrated': '#059669',
        }
        for name, (fpr, tpr, auc_val) in curves.items():
            color = palette.get(name, None)
            plt.plot(fpr, tpr, label=f'{name} (AUC={auc_val:.3f})', linewidth=2, color=color)
        plt.plot([0, 1], [0, 1], '--', color='gray', label='Chance')
        plt.xlim([0.0, 1.0])
        plt.ylim([0.0, 1.05])
        plt.xlabel('False Positive Rate (1 - Specificity)')
        plt.ylabel('True Positive Rate (Recall)')
        plt.title('ROC Comparison — Hypoxia Classification')
        plt.legend(loc='lower right')
        plt.grid(alpha=0.2)
        plt.tight_layout()
        plt.show()
    else:
        print('⚠️ No curves to plot; check that model files exist in models directory.')

# Summarize error patterns for the best model
from sklearn.metrics import confusion_matrix

# Use calibrated probabilities if available, otherwise raw
try:
    y_score_err = cal_proba
except NameError:
    if hasattr(best_model, 'predict_proba'):
        y_score_err = best_model.predict_proba(X_test)[:, 1]
    elif hasattr(best_model, 'decision_function'):
        y_score_err = best_model.decision_function(X_test)
    else:
        y_score_err = None

y_pred_err = best_model.predict(X_test)

cm = confusion_matrix(y_test, y_pred_err)
print('Confusion matrix (rows=true, cols=pred):\n', cm)

# Build a small analysis table
err_df = X_test_df.copy()
err_df['true'] = y_test.values
err_df['pred'] = y_pred_err
err_df['score'] = y_score_err if y_score_err is not None else np.nan

false_neg = err_df[(err_df['true'] == 1) & (err_df['pred'] == 0)]
false_pos = err_df[(err_df['true'] == 0) & (err_df['pred'] == 1)]
true_pos = err_df[(err_df['true'] == 1) & (err_df['pred'] == 1)]
true_neg = err_df[(err_df['true'] == 0) & (err_df['pred'] == 0)]

print(f'Counts — TP: {len(true_pos)}, FP: {len(false_pos)}, FN: {len(false_neg)}, TN: {len(true_neg)}')

# Compare feature means for FN vs TP (where we most care)
if len(false_neg) > 0 and len(true_pos) > 0:
    fn_tp_means = pd.DataFrame({
        'FN_mean': false_neg[feature_columns].mean(),
        'TP_mean': true_pos[feature_columns].mean(),
    }).sort_values('FN_mean', ascending=False)
    print('\nFeature means (FN vs TP) — watch for large gaps that may indicate miscalibration or missing interactions:')
    print(fn_tp_means.head(8))
else:
    print('No false negatives to summarize.')

# Quick glimpse of high-score FPs (if any)
if len(false_pos) > 0:
    top_fp = false_pos.sort_values('score', ascending=False).head(5)[feature_columns + ['score']]
    print('\nTop-scoring false positives (feature snapshot):')
    print(top_fp)
else:
    print('No false positives to inspect.')

Confusion matrix (rows=true, cols=pred):
 [[315  44]
 [ 17 342]]
Counts — TP: 342, FP: 44, FN: 17, TN: 315

Feature means (FN vs TP) — watch for large gaps that may indicate miscalibration or missing interactions:
            FN_mean     TP_mean
poc      527.646083  496.774446
sst       29.315045   29.798069
chlor_a   10.065308    9.160466
nflh       0.346490    0.259299
Rrs_547    0.009437    0.007357
Rrs_555    0.009366    0.007219
Rrs_531    0.008525    0.006824
Rrs_645    0.006665    0.003963

Top-scoring false positives (feature snapshot):
         chlor_a      nflh         poc        sst   Rrs_412   Rrs_443  \
99155   2.844585  0.286568  336.566650  30.149412  0.002204  0.002776   
99558   7.602374  0.164738  489.479981  29.863333  0.001752  0.002287   
30943  12.548010  0.135143  534.666687  31.475262  0.003642  0.004211   
31083  12.548010  0.135143  534.666687  31.475262  0.003642  0.004211   
30914  12.548010  0.135143  534.666687  31.475262  0.003642  0.004211   

        Rrs_469   Rrs_488   Rrs_531   Rrs_547   Rrs_555   Rrs_645   Rrs_667  \
99155  0.003416  0.003844  0.004506  0.004420  0.004229  0.001052  0.000934   
99558  0.002890  0.003317  0.004467  0.004649  0.004542  0.001565  0.001382   
30943  0.005196  0.005676  0.009309  0.010499  0.010354  0.007108  0.006409   
31083  0.005196  0.005676  0.009309  0.010499  0.010354  0.007108  0.006409   
30914  0.005196  0.005676  0.009309  0.010499  0.010354  0.007108  0.006409   

        Rrs_678     score  
99155  0.001045  0.942652  
99558  0.001445  0.937517  
30943  0.006021  0.937316  
31083  0.006021  0.937316  
30914  0.006021  0.937316

# Stratified k-fold CV metrics for robustness
from sklearn.model_selection import StratifiedKFold, cross_val_score

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

scoring = {
    'f1': 'f1',
    'balanced_accuracy': 'balanced_accuracy',
    'roc_auc': 'roc_auc',
}

cv_scores = {}
for metric_name, metric in scoring.items():
    scores = cross_val_score(best_model, X, y, cv=cv, scoring=metric, n_jobs=-1)
    cv_scores[metric_name] = scores
    print(f'{metric_name}: mean={scores.mean():.3f}, std={scores.std():.3f}, min={scores.min():.3f}, max={scores.max():.3f}')

# Quick table
cv_df = pd.DataFrame(cv_scores)
cv_df

f1: mean=0.919, std=0.006, min=0.912, max=0.928
balanced_accuracy: mean=0.916, std=0.006, min=0.908, max=0.925
roc_auc: mean=0.960, std=0.003, min=0.955, max=0.963

# Apply trained model to satellite raster files for spatial prediction
# Input raster files (MODIS GeoTIFFs with 14 bands)
# files D:\my works\bloom\data\MODIS
rasters=glob.glob(r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\modis out\*.tif')
# path
output_path=r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
os.makedirs(output_path, exist_ok=True)

# Load the fitted scaler used during model training
scaler_path = os.path.join(models_directory, 'MinMaxScaler.pkl')
scaler = joblib.load(scaler_path)
print(f'✓ Loaded fitted scaler from: {scaler_path}')

print(f'Found {len(rasters)} raster files to process.')

# Process each raster
for raster in rasters:
    print(f'\nProcessing: {os.path.basename(raster)}')
    
    with rasterio.open(raster) as src:
        # Read all bands into a numpy array (bands, height, width)
        array = src.read()
        # Capture spatial metadata for output
        meta = src.meta
    
    print(f'  Raster shape (bands, height, width): {array.shape}')
    
    # Stack bands into (height, width, bands) format
    array_stack = np.stack(array, axis=-1)
    print(f'  Stacked shape (height, width, bands): {array_stack.shape}')
    
    # Reshape to 2D for prediction: (n_pixels, n_bands)
    array_reshaped = array_stack.reshape(-1, array_stack.shape[-1])
    print(f'  Reshaped for model (n_pixels, n_bands): {array_reshaped.shape}')
    
    # Create DataFrame with proper feature names
    df_array = pd.DataFrame(
        array_reshaped,
        columns=['chlor_a', 'nflh', 'poc', 'sst', 'Rrs_412',
                 'Rrs_443', 'Rrs_469', 'Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555',
                 'Rrs_645', 'Rrs_667', 'Rrs_678']
    )
    
    print(f'  DataFrame shape: {df_array.shape}')
    print(f'  First 2 rows:\n{df_array.head(2)}')
    
    # Normalize features using the same scaler used during training
    df_array_normalized = scaler.transform(df_array)
    print(f'  Normalized shape: {df_array_normalized.shape}')
    
    # Load the best trained model
    best_model_path = os.path.join(models_directory, 'best_model.pkl')
    model = joblib.load(best_model_path)
    print(f'  Loaded model from: {best_model_path}')
    
    # Generate predictions for all pixels
    predictions = model.predict(df_array_normalized)
    print(f'  Predictions shape: {predictions.shape}')
    
    # Reshape predictions back to 2D raster (height, width)
    predictions_2d = predictions.reshape(array_stack.shape[0], array_stack.shape[1])
    print(f'  Predictions 2D shape (height, width): {predictions_2d.shape}')
    
    # Update metadata for single-band output raster
    meta.update(count=1, dtype=rasterio.uint8, height=predictions_2d.shape[0], width=predictions_2d.shape[1])
    
    # Save prediction as GeoTIFF with spatial metadata
    output_filename = os.path.splitext(os.path.basename(raster))[0] + '_hypoxia_pred.tif'
    output_file = os.path.join(output_path, output_filename)
    
    with rasterio.open(output_file, 'w', **meta) as dst:
        dst.write(predictions_2d.astype(rasterio.uint8), 1)
    
    print(f'  ✓ Saved prediction to: {output_file}')
    
    # Clear memory
    del array, array_stack, array_reshaped, df_array, df_array_normalized, predictions, predictions_2d

print(f'\n✓ Raster prediction complete! All outputs saved to: {output_path}')

import re

hurrican = pd.read_csv(r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Hurricane.csv')
# Convert hurricane Date strings to datetime (supports formats like 'May-2010' or '01-May-2010')
def convert_date(date):
    for fmt in ('%b-%Y', '%d-%b-%Y'):
        try:
            return pd.to_datetime(date, format=fmt)
        except ValueError:
            continue
    return pd.NaT  # Return not-a-time for invalid dates

hurrican['Date'] = hurrican['Date'].apply(convert_date)

path = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
# read csv
df_area = pd.read_csv(os.path.join(path, 'area.csv'))

# Helper to parse file names like 'SMI2010_05' into timestamps (defaults to day=1)
def parse_file_date(val: str):
    # First, try direct parsing (handles ISO/mixed)
    dt = pd.to_datetime(val, errors='coerce', format=None)
    if pd.notnull(dt):
        return dt
    # Fallback: extract YYYY and MM from patterns like SMI2010_05 or 2010-05
    m = re.search(r'(\d{4})[_-]?(\d{2})', str(val))
    if m:
        year = int(m.group(1))
        month = int(m.group(2))
        return pd.Timestamp(year=year, month=month, day=1)
    return pd.NaT

# Select hypoxia class and parse dates
df_hypo = df_area[['file', 'class_1_area_sqkm']].copy()
df_hypo['file_dt'] = df_hypo['file'].apply(parse_file_date)
# Drop rows where date could not be parsed
df_hypo = df_hypo.dropna(subset=['file_dt'])

# Sort by parsed date
df_hypo.sort_values(by='file_dt', inplace=True)

# Temporal features
df_hypo['month'] = df_hypo['file_dt'].dt.month
df_hypo['season'] = df_hypo['month'] % 12 + 1
df_hypo['season'] = df_hypo['season'].apply(
    lambda x: 'Spring' if 3 <= x <= 5 else (
        'Summer' if 6 <= x <= 8 else (
            'Winter' if x == 12 or x <= 2 else None)))

# Remove rows with season as None (Autumn months)
df_hypo = df_hypo[df_hypo['season'].notna()]
palette = sns.color_palette("hsv", 12)

# Create a figure with three subplots
fig, (ax3, ax2, ax1) = plt.subplots(3, 1, figsize=(15, 15))

# Monthly boxplot
sns.boxplot(x='month', y='class_1_area_sqkm', data=df_hypo, ax=ax3, palette=palette)
ax3.set_ylabel('Predicted Hypoxic area over sqkm')
ax3.set_xlabel('Month')
ax3.annotate('(a)', xy=(0.01, 0.95), xycoords='axes fraction', fontsize=12,
             horizontalalignment='left', verticalalignment='top')

# Seasonal violin + swarm
sns.violinplot(x='season', y='class_1_area_sqkm', data=df_hypo, ax=ax2, palette='Set3', inner=None)
sns.swarmplot(x='season', y='class_1_area_sqkm', data=df_hypo, ax=ax2, color='k', size=3)
ax2.set_ylabel('Predicted Hypoxic area over sqkm')
ax2.annotate('(b)', xy=(0.01, 0.95), xycoords='axes fraction', fontsize=12,
             horizontalalignment='left', verticalalignment='top')

# Time series with trend and optional hurricane markers
sns.lineplot(x=df_hypo['file_dt'], y=df_hypo['class_1_area_sqkm'], ax=ax1, color='blue')
if 'class_1_area_sqkm' in hurrican.columns:
    sns.scatterplot(x=hurrican['Date'], y=hurrican['class_1_area_sqkm'], color='red', s=100,
                    ax=ax1, label='Hurricane Events')
ax1.set_xlim(df_hypo['file_dt'].min(), df_hypo['file_dt'].max())
z = np.polyfit(range(len(df_hypo)), df_hypo['class_1_area_sqkm'], 1)
p = np.poly1d(z)
ax1.plot(df_hypo['file_dt'], p(range(len(df_hypo))), "r--", label='Trend line')
mean_val = df_hypo['class_1_area_sqkm'].mean()
ax1.axhline(mean_val, color='gray', linestyle='dotted', label='Averaged line')
ax1.set_ylabel('Predicted Hypoxic area over sqkm')
ax1.set_xlabel('')
ax1.annotate('(c)', xy=(0.01, 0.95), xycoords='axes fraction', fontsize=12,
             horizontalalignment='left', verticalalignment='top')
ax1.legend(loc='upper center')

plt.tight_layout()
# plt.savefig(r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\figure\seasonal_analysis_v3.jpg', dpi=600, bbox_inches='tight')
plt.show()

C:\Users\hafez\AppData\Local\Temp\ipykernel_59372\4094274612.py:58: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  sns.boxplot(x='month', y='class_1_area_sqkm', data=df_hypo, ax=ax3, palette=palette)
C:\Users\hafez\AppData\Local\Temp\ipykernel_59372\4094274612.py:58: UserWarning: The palette list has more values (12) than needed (9), which may not be intended.
  sns.boxplot(x='month', y='class_1_area_sqkm', data=df_hypo, ax=ax3, palette=palette)
C:\Users\hafez\AppData\Local\Temp\ipykernel_59372\4094274612.py:65: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  sns.violinplot(x='season', y='class_1_area_sqkm', data=df_hypo, ax=ax2, palette='Set3', inner=None)

# Get probability of hypoxia (class 1) for each pixel
probabilities = model.predict_proba(df_array_normalized)[:, 1]

# Generate uncertainty (probability) rasters alongside binary predictions
# Input raster files (MODIS GeoTIFFs with 14 bands)
rasters = glob.glob(r'C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\modis out\*.tif')

# Output directories
output_path_binary = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\probability"
output_path_uncertainty = os.path.join(output_path_binary, "uncertainty")
os.makedirs(output_path_binary, exist_ok=True)
os.makedirs(output_path_uncertainty, exist_ok=True)

print(f'Found {len(rasters)} raster files to process for uncertainty quantification.')

# Process each raster
for raster in rasters:
    print(f'\nProcessing: {os.path.basename(raster)}')
    
    with rasterio.open(raster) as src:
        # Read all bands into a numpy array (bands, height, width)
        array = src.read()
        # Capture spatial metadata for output
        meta = src.meta
    
    print(f'  Raster shape (bands, height, width): {array.shape}')
    
    # Stack bands into (height, width, bands) format
    array_stack = np.stack(array, axis=-1)
    
    # Reshape to 2D for prediction: (n_pixels, n_bands)
    array_reshaped = array_stack.reshape(-1, array_stack.shape[-1])
    
    # Create DataFrame with proper feature names
    df_array = pd.DataFrame(
        array_reshaped,
        columns=['chlor_a', 'nflh', 'poc', 'sst', 'Rrs_412',
                 'Rrs_443', 'Rrs_469', 'Rrs_488', 'Rrs_531', 'Rrs_547', 'Rrs_555',
                 'Rrs_645', 'Rrs_667', 'Rrs_678']
    )
    
    # Normalize features using the same scaler used during training
    df_array_normalized = scaler.transform(df_array)
    
    # Load the best trained model
    best_model_path = os.path.join(models_directory, 'best_model.pkl')
    model = joblib.load(best_model_path)
    
    # Generate binary predictions
    predictions = model.predict(df_array_normalized)
    
    # Generate probability predictions (use calibrated model if available)
    try:
        # Try to load and use the calibrated model for better uncertainty estimates
        calibrated_model_path = os.path.join(models_directory, f'{best_model_name}_calibrated.pkl')
        if os.path.exists(calibrated_model_path):
            calibrator_loaded = joblib.load(calibrated_model_path)
            probabilities = calibrator_loaded.predict_proba(df_array_normalized)[:, 1]
            print(f'  ✓ Using calibrated probabilities for uncertainty')
        else:
            # Fall back to raw probabilities from best model
            if hasattr(model, 'predict_proba'):
                probabilities = model.predict_proba(df_array_normalized)[:, 1]
            else:
                probabilities = model.decision_function(df_array_normalized)
                # Normalize decision function to [0, 1] range
                probabilities = (probabilities - probabilities.min()) / (probabilities.max() - probabilities.min() + 1e-9)
            print(f'  Using raw model probabilities for uncertainty')
    except Exception as e:
        print(f'  Warning: Could not load calibrated model, using raw probabilities: {e}')
        if hasattr(model, 'predict_proba'):
            probabilities = model.predict_proba(df_array_normalized)[:, 1]
        else:
            probabilities = model.decision_function(df_array_normalized)
            probabilities = (probabilities - probabilities.min()) / (probabilities.max() - probabilities.min() + 1e-9)
    
    # Reshape predictions and probabilities back to 2D raster (height, width)
    predictions_2d = predictions.reshape(array_stack.shape[0], array_stack.shape[1])
    probabilities_2d = probabilities.reshape(array_stack.shape[0], array_stack.shape[1])
    
    print(f'  Predictions 2D shape: {predictions_2d.shape}')
    print(f'  Probabilities range: [{probabilities_2d.min():.3f}, {probabilities_2d.max():.3f}]')
    
    # Save binary prediction raster
    meta_binary = meta.copy()
    meta_binary.update(count=1, dtype=rasterio.uint8, height=predictions_2d.shape[0], width=predictions_2d.shape[1])
    
    output_filename_binary = os.path.splitext(os.path.basename(raster))[0] + '_hypoxia_pred.tif'
    output_file_binary = os.path.join(output_path_binary, output_filename_binary)
    
    with rasterio.open(output_file_binary, 'w', **meta_binary) as dst:
        dst.write(predictions_2d.astype(rasterio.uint8), 1)
    print(f'  ✓ Saved binary prediction to: {output_file_binary}')
    
    # Save uncertainty (probability) raster as uint8 (0-255 scale)
    # 0 = certain non-hypoxic, 255 = certain hypoxic, 127 = maximum uncertainty
    probabilities_uint8 = (probabilities_2d * 255).astype(rasterio.uint8)
    
    meta_uncertainty = meta.copy()
    meta_uncertainty.update(count=1, dtype=rasterio.uint8, height=probabilities_2d.shape[0], width=probabilities_2d.shape[1])
    
    output_filename_uncertainty = os.path.splitext(os.path.basename(raster))[0] + '_hypoxia_uncertainty.tif'
    output_file_uncertainty = os.path.join(output_path_uncertainty, output_filename_uncertainty)
    
    with rasterio.open(output_file_uncertainty, 'w', **meta_uncertainty) as dst:
        dst.write(probabilities_uint8, 1)
    print(f'  ✓ Saved uncertainty raster to: {output_file_uncertainty}')
    
    # Compute and display uncertainty stats
    entropy = -probabilities_2d * np.log2(probabilities_2d + 1e-9) - (1 - probabilities_2d) * np.log2(1 - probabilities_2d + 1e-9)
    high_confidence = (probabilities_2d > 0.8) | (probabilities_2d < 0.2)
    pct_confident = 100 * high_confidence.sum() / high_confidence.size
    print(f'  High-confidence pixels (prob > 0.8 or < 0.2): {pct_confident:.1f}%')
    
    # Clear memory
    del array, array_stack, array_reshaped, df_array, df_array_normalized, predictions, probabilities
    del predictions_2d, probabilities_2d, probabilities_uint8, entropy, high_confidence

print(f'\n✓ Uncertainty quantification complete!')
print(f'  Binary predictions saved to: {output_path_binary}')
print(f'  Uncertainty rasters saved to: {output_path_uncertainty}')

# Visualize a single uncertainty raster (uint8 0–255 → probability 0–1)
import os
import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import plotting_extent

uncertainty_path = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\probability\uncertainty\OCEANDATA_MODIS-Aqua_L3SMI2002-07-01_hypoxia_uncertainty.tif"

with rasterio.open(uncertainty_path) as src:
    data_u8 = src.read(1)  # uint8 values 0–255
    extent = plotting_extent(src)
    res = src.res
    crs = src.crs
    bounds = src.bounds

# Convert to probability scale (0–1)
data_prob = data_u8.astype(np.float32) / 255.0
uncertainty = 1.0 - np.minimum(np.abs(data_prob - 0.5) * 2.0, 1.0)
# Plot probability map and histogram

try:
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    HAS_CARTOPY = True
except Exception:
    HAS_CARTOPY = False


# Optional cartopy overlay (coastline)
if HAS_CARTOPY:
    fig, ax = plt.subplots(1, 1, figsize=(10, 6), subplot_kw={'projection': ccrs.PlateCarree()})
    im = ax.imshow(uncertainty, extent=extent, cmap='YlOrBr', vmin=0, vmax=1, transform=ccrs.PlateCarree())
    ax.add_feature(cfeature.LAND, facecolor='lightgrey', zorder=0)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.8, edgecolor='black', zorder=2)
    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False
    plt.colorbar(im, ax=ax, orientation='vertical', label='Uncertainty (0–1)')
    ax.set_title('Uncertainty (with coastlines)', fontsize=12)
    plt.tight_layout()
    plt.show()

from rasterio.warp import calculate_default_transform, reproject, Resampling


path=r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
#find .tif
rasters=glob.glob(os.path.join(path,'*.tif'))

# read first
with rasterio.open(rasters[0]) as src:
    array = src.read(1)
    meta = src.meta
    resolution = src.res
    # crs
    crs = src.crs

    # Define the target CRS
    dst_crs = 'EPSG:3814'

    # Calculate the transform and the new dimensions
    transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)

    # Create a new array to hold the reprojected data
    dst_array = np.empty((height, width))

    # Reproject the data
    reproject(
        source=array,
        destination=dst_array,
        src_transform=src.transform,
        src_crs=src.crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest)

# print array and resolution
print(f"Original array shape: {array.shape}")
print(f"Reprojected array shape: {dst_array.shape}")
print(f"Original Resolution: {resolution}")
print(f"Original CRS: {crs}")
print(f"New CRS: {dst_crs}")
# reprojected resolution
dst_resolution = (transform.a, transform.e)
print(f"Reprojected Resolution: {dst_resolution}")
# reprojected resolution
dst_resolution = (transform.a, transform.e)
print(f"Reprojected Resolution: {dst_resolution}")
# Assuming 0 and 1 are the classes in your raster
# Assuming 0 and 1 are the classes in your raster
class_0_area_sqm = np.count_nonzero(dst_array == 0) * (dst_resolution[0] * -dst_resolution[1])
class_1_area_sqm = np.count_nonzero(dst_array == 1) * (dst_resolution[0] * -dst_resolution[1])

# Convert to square kilometers
class_0_area_sqkm = class_0_area_sqm / 1e6
class_1_area_sqkm = class_1_area_sqm / 1e6

print(f"Area of class 0: {class_0_area_sqkm} square kilometers")
print(f"Area of class 1: {class_1_area_sqkm} square kilometers")
# calculate area of each file in rasters  and make dataframe: file name, class 0 area, class 1 area
df_area=[]
for raster in rasters:
    with rasterio.open(raster) as src:
        array = src.read(1)
        meta = src.meta
        resolution = src.res
        crs = src.crs
        dst_crs = 'EPSG:3814'
        transform, width, height = calculate_default_transform(src.crs, dst_crs, src.width, src.height, *src.bounds)
        dst_array = np.empty((height, width))
        reproject(
            source=array,
            destination=dst_array,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=transform,
            dst_crs=dst_crs,
            resampling=Resampling.nearest)
        dst_resolution = (transform.a, transform.e)
        class_0_area_sqm = np.count_nonzero(dst_array == 0) * (dst_resolution[0] * -dst_resolution[1])
        class_1_area_sqm = np.count_nonzero(dst_array == 1) * (dst_resolution[0] * -dst_resolution[1])
        class_0_area_sqkm = class_0_area_sqm / 1e6
        class_1_area_sqkm = class_1_area_sqm / 1e6
        name=raster[-18:-8]
        # make dataframe, then append to list
        df_temp=pd.DataFrame({'file':[name],'class_0_area_sqkm':[class_0_area_sqkm],'class_1_area_sqkm':[class_1_area_sqkm]})
        df_area.append(df_temp)
        del df_temp, array, dst_array, class_0_area_sqm, class_1_area_sqm, class_0_area_sqkm, class_1_area_sqkm
# convert list to dataframe
df_area=pd.concat(df_area,ignore_index=True)
# save to csv
df_area.to_csv(os.path.join(path,'area.csv'),index=False)
df_area

Original array shape: (850, 1936)
Reprojected array shape: (976, 1960)
Original Resolution: (0.008983152841195214, 0.008983152841195214)
Original CRS: EPSG:4326
New CRS: EPSG:3814
Reprojected Resolution: (906.05699435686, -906.05699435686)
Reprojected Resolution: (906.05699435686, -906.05699435686)
Area of class 0: 1461363.0373606663 square kilometers
Area of class 1: 109060.96201322679 square kilometers

import matplotlib.pyplot as plt
import rasterio
import rasterio.plot
import numpy as np
import glob
import os

path = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
saved_file=r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\averagedPredicted"
fig, axs = plt.subplots(5, 2, figsize=(15, 30), constrained_layout=True)  # Create a figure with 10 subplots
axs = axs.ravel()  # Flatten the array of axes

# Define a variable to store the maximum intensity of the hypoxic zone across all plots
max_intensity = 0

for i, year in enumerate(range(2011, 2021)):
    # Select only files containing the current year
    rasters = glob.glob(os.path.join(path, f'*{year}*.tif'))

    # Read all files and create an array stack
    array_stack = []
    for raster in rasters:
        with rasterio.open(raster) as src:
            array = src.read(1)
            array_stack.append(array)

    # Stack arrays along the first axis
    array_stack = np.stack(array_stack, axis=0)

    # Calculate the average
    array_avg = np.nanmean(array_stack, axis=0)  # Use np.nanmean to ignore NaN values
    array_avg[array_avg == 0] = np.nan  # Remove 0 values

    # Extract geospatial information
    with rasterio.open(rasters[0]) as src:
        transform = src.transform
        crs = src.crs

    # Plot on the current subplot
    im = axs[i].imshow(array_avg, extent=rasterio.plot.plotting_extent(src), cmap='viridis')
    axs[i].set_xlabel('Longitude')
    axs[i].set_ylabel('Latitude')
    axs[i].set_title(f'Average of predicted hypoxic area ({year})', fontsize=14, fontweight='bold')

    # Update max_intensity if necessary
    max_intensity = max(max_intensity, np.nanmax(array_avg))

# Add a colorbar to the figure
cbar = fig.colorbar(im, ax=axs, orientation='horizontal', pad=0.02)
cbar.set_label('Intensity of Hypoxic Zone')
# save fig
#plt.savefig(r'D:\my works\bloom\figure\average_raster_data.jpg',dpi=600,bbox_inches='tight')
plt.show()
# clean memory
del array_stack, array_avg, transform, crs, im, cbar

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from mpl_toolkits.axes_grid1 import make_axes_locatable

path = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
saved_file = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\averagedPredicted"
os.makedirs(saved_file, exist_ok=True)

fig, axs = plt.subplots(4, 3, figsize=(18, 12),
                        gridspec_kw={'wspace': 0.05, 'hspace': 0.05},
                        subplot_kw={'projection': ccrs.PlateCarree()})
axs = axs.ravel()

years = list(range(2011, 2021))[:12]  # 12 years for 3x4 grid

im = None  # Initialize im for colorbar

for i, year in enumerate(years):
    rasters = glob.glob(os.path.join(path, f'*{year}*.tif'))
    if not rasters:
        axs[i].axis('off')
        continue

    array_stack = []
    for raster in rasters:
        with rasterio.open(raster) as src:
            array = src.read(1)
            array_stack.append(array)

    array_stack = np.stack(array_stack, axis=0)
    array_avg = np.nanmean(array_stack, axis=0)
    array_avg[array_avg == 0] = np.nan

    with rasterio.open(rasters[0]) as src:
        extent = rasterio.plot.plotting_extent(src)
        im = axs[i].imshow(array_avg, extent=extent, cmap='jet', transform=ccrs.PlateCarree())
        #axs[i].add_feature(cfeature.COASTLINE, linewidth=0.8, edgecolor='black', zorder=2)
        axs[i].add_feature(cfeature.LAND, facecolor='lightgrey', zorder=0)
        axs[i].annotate(f'({chr(97+i)})', xy=(0.01, 0.97), xycoords='axes fraction',
                        fontsize=14, fontweight='bold', ha='left', va='top')

        # Gridlines: only left and bottom, no right/top labels
        gl = axs[i].gridlines(draw_labels=True, linewidth=0.5, color='gray',
                              alpha=0.5, linestyle='--', zorder=3)
        gl.top_labels = False
        gl.right_labels = False
        gl.xlabel_style = {'size': 10}
        gl.ylabel_style = {'size': 10}
        gl.left_labels = (i % 3 == 0)
        gl.bottom_labels = (i // 3 == 3)
        if not gl.left_labels:
            axs[i].set_yticks([])
        if not gl.bottom_labels:
            axs[i].set_xticks([])

# Hide any unused subplots
for j in range(len(years), len(axs)):
    axs[j].axis('off')

# Place colorbar in the last two empty subplots (axs[9] and axs[10])
from matplotlib.gridspec import GridSpec
gs = axs[0].get_gridspec()
# Remove the last two axes to make space for colorbar


# Create a new axis spanning the last row, columns 1 and 2
cax = fig.add_subplot(gs[3, 1:3])
cbar = fig.colorbar(im, cax=cax, orientation='horizontal', aspect=30)
cbar.set_label('Intensity of Hypoxic Zone', fontsize=14, fontweight='bold')
cax.set_aspect(0.08)  # Lower value = thinner
# save fig
plt.savefig(os.path.join(saved_file, 'average_predicted_hypoxic_2011_2020.jpg'), dpi=600, bbox_inches='tight')
plt.show()
del array_stack, array_avg, im, cbar
# caption: Average predicted hypoxic zone intensity from 2011 to 2020:Here

# Example: Simple Dashboard Prototype with Plotly Dash
# Install: pip install dash plotly rasterio

import dash
from dash import dcc, html, Input, Output
import plotly.graph_objects as go
import rasterio
import numpy as np
import glob
import os

# Initialize Dash app
app = dash.Dash(__name__)

# Load prediction rasters
raster_dir = r"C:\Users\hafez\MSU\Research\msGOM\mssound\bloom\data\Predicted rast"
pred_files = sorted(glob.glob(os.path.join(raster_dir, "*_hypoxia_pred.tif")))
uncertainty_files = sorted(glob.glob(os.path.join(raster_dir, "*_hypoxia_uncertainty.tif")))

# Extract dates from filenames for dropdown
dates = [os.path.basename(f)[-18:-8] for f in pred_files]

# Dashboard layout
app.layout = html.Div([
    html.H1("Gulf of Mexico Hypoxia Prediction Dashboard", style={'textAlign': 'center'}),
    
    html.Div([
        html.Label("Select Date:"),
        dcc.Dropdown(
            id='date-dropdown',
            options=[{'label': date, 'value': i} for i, date in enumerate(dates)],
            value=0,
            style={'width': '300px'}
        ),
    ], style={'padding': '20px'}),
    
    html.Div([
        html.Div([
            html.H3("Hypoxia Prediction (Binary)"),
            dcc.Graph(id='prediction-map')
        ], style={'width': '48%', 'display': 'inline-block'}),
        
        html.Div([
            html.H3("Prediction Uncertainty (Confidence)"),
            dcc.Graph(id='uncertainty-map')
        ], style={'width': '48%', 'display': 'inline-block', 'float': 'right'})
    ]),
    
    html.Div([
        html.H3("Statistics"),
        html.Div(id='stats-panel', style={'fontSize': '16px', 'padding': '20px'})
    ])
])

# Callback to update maps and statistics
@app.callback(
    [Output('prediction-map', 'figure'),
     Output('uncertainty-map', 'figure'),
     Output('stats-panel', 'children')],
    [Input('date-dropdown', 'value')]
)
def update_dashboard(selected_idx):
    # Load prediction raster
    with rasterio.open(pred_files[selected_idx]) as src:
        pred_data = src.read(1)
        bounds = src.bounds
    
    # Load uncertainty raster
    with rasterio.open(uncertainty_files[selected_idx]) as src:
        uncert_data = src.read(1).astype(float)
    
    # Create prediction heatmap
    pred_fig = go.Figure(data=go.Heatmap(
        z=pred_data,
        colorscale=[[0, 'lightblue'], [1, 'red']],
        colorbar=dict(title="Class", tickvals=[0, 1], ticktext=["Non-hypoxic", "Hypoxic"])
    ))
    pred_fig.update_layout(
        xaxis_title="Longitude",
        yaxis_title="Latitude",
        height=500
    )
    
    # Create uncertainty heatmap
    uncert_fig = go.Figure(data=go.Heatmap(
        z=uncert_data,
        colorscale='Viridis',
        colorbar=dict(title="Confidence (0-255)")
    ))
    uncert_fig.update_layout(
        xaxis_title="Longitude",
        yaxis_title="Latitude",
        height=500
    )
    
    # Compute statistics
    total_pixels = np.sum(~np.isnan(pred_data))
    hypoxic_pixels = np.sum(pred_data == 1)
    hypoxic_pct = (hypoxic_pixels / total_pixels * 100) if total_pixels > 0 else 0
    avg_uncertainty = np.nanmean(uncert_data)
    high_conf_pixels = np.sum((uncert_data > 200) | (uncert_data < 55))
    high_conf_pct = (high_conf_pixels / total_pixels * 100) if total_pixels > 0 else 0
    
    stats_text = [
        html.P(f"📅 Date: {dates[selected_idx]}"),
        html.P(f"🌊 Hypoxic Area: {hypoxic_pct:.1f}% of pixels ({hypoxic_pixels:,} pixels)"),
        html.P(f"📊 Average Uncertainty: {avg_uncertainty:.1f} (0=low, 255=high)"),
        html.P(f"✅ High Confidence Pixels: {high_conf_pct:.1f}% (>{200} or <{55})")
    ]
    
    return pred_fig, uncert_fig, stats_text

# Note: To run dashboard, execute: app.run_server(debug=True, port=8050)
# Then open browser to http://localhost:8050
print("Dashboard code ready. To launch, uncomment: app.run_server(debug=True, port=8050)")

# Example: FastAPI Model Serving Endpoint
# Save as api.py and run: uvicorn api:app --reload
# Install: pip install fastapi uvicorn pydantic joblib scikit-learn

"""
from fastapi import FastAPI, HTTPException
from pydantic import BaseModel
import joblib
import numpy as np
import pandas as pd
from typing import List

app = FastAPI(title="Hypoxia Prediction API", version="1.0")

# Load model and scaler at startup
model = joblib.load('models/best_model.pkl')
scaler = joblib.load('models/scaler.pkl')

# Define input schema
class PredictionInput(BaseModel):
    chlor_a: float
    nflh: float
    poc: float
    sst: float
    Rrs_412: float
    Rrs_443: float
    Rrs_469: float
    Rrs_488: float
    Rrs_531: float
    Rrs_547: float
    Rrs_555: float
    Rrs_645: float
    Rrs_667: float
    Rrs_678: float

class BatchPredictionInput(BaseModel):
    features: List[PredictionInput]

# Health check endpoint
@app.get("/")
def health_check():
    return {"status": "healthy", "model": "best_model.pkl"}

# Single prediction endpoint
@app.post("/predict")
def predict(input_data: PredictionInput):
    try:
        # Convert input to DataFrame
        features_df = pd.DataFrame([input_data.dict()])
        
        # Scale features
        features_scaled = scaler.transform(features_df)
        
        # Make prediction
        prediction = model.predict(features_scaled)[0]
        probability = model.predict_proba(features_scaled)[0, 1]
        
        return {
            "prediction": int(prediction),
            "probability_hypoxic": float(probability),
            "class_label": "Hypoxic" if prediction == 1 else "Non-hypoxic"
        }
    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))

# Batch prediction endpoint
@app.post("/predict_batch")
def predict_batch(batch_input: BatchPredictionInput):
    try:
        # Convert batch to DataFrame
        features_df = pd.DataFrame([item.dict() for item in batch_input.features])
        
        # Scale features
        features_scaled = scaler.transform(features_df)
        
        # Make predictions
        predictions = model.predict(features_scaled)
        probabilities = model.predict_proba(features_scaled)[:, 1]
        
        results = [
            {
                "prediction": int(pred),
                "probability_hypoxic": float(prob),
                "class_label": "Hypoxic" if pred == 1 else "Non-hypoxic"
            }
            for pred, prob in zip(predictions, probabilities)
        ]
        
        return {"predictions": results, "count": len(results)}
    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))

# Model metadata endpoint
@app.get("/model_info")
def model_info():
    return {
        "model_type": str(type(model).__name__),
        "input_features": 14,
        "feature_names": [
            "chlor_a", "nflh", "poc", "sst",
            "Rrs_412", "Rrs_443", "Rrs_469", "Rrs_488",
            "Rrs_531", "Rrs_547", "Rrs_555", "Rrs_645",
            "Rrs_667", "Rrs_678"
        ],
        "output_classes": ["Non-hypoxic", "Hypoxic"]
    }
"""

print("FastAPI code ready. Save to api.py and run: uvicorn api:app --reload")
print("API will be available at: http://localhost:8000")
print("Interactive docs: http://localhost:8000/docs")

FROM python:3.11-slim
WORKDIR /app
COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt
COPY models/ ./models/
COPY api.py .
EXPOSE 8000
CMD ["uvicorn", "api:app", "--host", "0.0.0.0", "--port", "8000"]

	latitude	longitude	activitystartdate	resultmeasurevalue_mg_l
0	30.2796	-87.5417	1999-10-04 00:00:00+00:00	6.3
1	30.2796	-87.5417	1999-11-01 00:00:00+00:00	6.8
2	30.2796	-87.5417	1999-12-06 00:00:00+00:00	7.4

	latitude	longitude	activitystartdate	resultmeasurevalue_mg_l	hypoxia	year
0	30.2796	-87.5417	1999-10-04 00:00:00+00:00	6.3	1	1999
1	30.2796	-87.5417	1999-11-01 00:00:00+00:00	6.8	1	1999
2	30.2796	-87.5417	1999-12-06 00:00:00+00:00	7.4	1	1999

	latitude	longitude	activitystartdate	resultmeasurevalue_mg_l	hypoxia	year	5_year_period	id
0	30.2796	-87.5417	1999-10-04 00:00:00+00:00	6.30	1	1999	1990	1
1	30.2796	-87.5417	1999-11-01 00:00:00+00:00	6.80	1	1999	1990	2
2	30.2796	-87.5417	1999-12-06 00:00:00+00:00	7.40	1	1999	1990	3
3	30.2796	-87.5417	2000-01-10 00:00:00+00:00	7.52	1	2000	2000	4
4	30.2796	-87.5417	2000-02-07 00:00:00+00:00	8.30	1	2000	2000	5
...	...	...	...	...	...	...	...	...
6725948	30.3808	-87.8320	2014-01-01 00:00:00+00:00	3.80	1	2014	2010	6093331
6725949	30.3808	-87.8320	2014-01-01 00:00:00+00:00	3.80	1	2014	2010	6093332
6725950	30.3808	-87.8320	2014-01-01 00:00:00+00:00	3.80	1	2014	2010	6093333
6725951	30.3808	-87.8320	2014-01-01 00:00:00+00:00	3.70	1	2014	2010	6093334
6725952	30.3808	-87.8320	2014-01-01 00:00:00+00:00	3.70	1	2014	2010	6093335

	resultmeasurevalue_mg_l	chlor_a	nflh	poc	sst	Rrs_412	Rrs_443	Rrs_469	Rrs_488	Rrs_531	Rrs_547	Rrs_555	Rrs_645	Rrs_667	Rrs_678
0	4.997	0.209567	0.019738	63.439999	25.814131	0.005095	0.005134	0.004986	0.004535	0.002386	0.001788	0.001522	0.000124	0.000172	0.000168
1	5.003	0.209567	0.019738	63.439999	25.814131	0.005095	0.005134	0.004986	0.004535	0.002386	0.001788	0.001522	0.000124	0.000172	0.000168

Model	Why Useful for Hypoxia	Advantage	Challenge
Deep Neural Networks (DNNs)	High-dimensional satellite & environmental data	Learn hierarchical non-linear features; handle temporal sequences	Requires large datasets (~10k+ samples); harder to interpret decisions
Graph Neural Networks (GNNs)	Spatial adjacency of hypoxic zones matters	Explicit representation of spatial correlations between ocean pixels	Complex implementation; requires spatial graph construction
Recurrent Neural Networks (RNNs/LSTMs)	Hypoxia has temporal memory (stratification persists)	Capture temporal dependencies (e.g., "SST trend + density stratification → risk")	Needs long time series per location; computationally expensive
Ensemble Boosting (XGBoost, LightGBM)	Handles mixed feature types & nonlinearity	Often outperforms classical RF in practice; built-in feature selection	Black-box nature; harder to interpret interactions
Uncertainty Quantification (Bayesian NNs, Variational Inference)	Coastal managers need prediction confidence bounds	Quantify uncertainty in hypoxia predictions; distinguish model vs. data uncertainty	High computational cost; requires posterior sampling
Transfer Learning (Pre-trained Vision Models)	Satellite imagery as spatial-spectral "images"	Leverage learned representations from MODIS/Sentinel time series; fine-tune on hypoxia	Domain-specific training data needed; generic vision models may not transfer well

	Classifier	Accuracy	Balanced_Accuracy	Precision	Recall	F1	ROC_AUC	Model_Path
0	Random Forest	0.915042	0.915042	0.886010	0.952646	0.918121	0.961946	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
6	XGBoost	0.913649	0.913649	0.881748	0.955432	0.917112	0.956755	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
3	K-Nearest Neighbors	0.909471	0.909471	0.880829	0.947075	0.912752	0.943603	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
4	Decision Tree	0.903900	0.903900	0.873711	0.944290	0.907631	0.940542	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
2	Support Vector Machine	0.884401	0.884401	0.828571	0.969359	0.893453	0.916570	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
7	MLP	0.881616	0.881616	0.834146	0.952646	0.889467	0.941167	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
1	Logistic Regression	0.876045	0.876045	0.806818	0.988858	0.888611	0.905638	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...
5	Naive Bayes	0.818942	0.818942	0.815427	0.824513	0.819945	0.874939	C:\Users\hafez\MSU\Research\msGOM\mssound\bloo...

	f1	balanced_accuracy	roc_auc
0	0.920892	0.918523	0.960992
1	0.928000	0.924939	0.963302
2	0.911647	0.907950	0.955025
3	0.914513	0.910042	0.961739
4	0.919028	0.916318	0.961415

Probability $p$	Entropy $H(p)$	Interpretation
$p = 0.0$ or $1.0$	$H = 0$	Maximum confidence (certain non-hypoxic or hypoxic)
$p = 0.5$	$H = 1$	Maximum uncertainty (equiprobable classes, decision boundary)
$p = 0.8$ or $0.2$	$H \approx 0.72$	Moderate confidence

Pixel Value (uint8)	Probability $p$	Meaning
0	0.00	Certain non-hypoxic
64	0.25	Low probability hypoxic
127	0.50	Maximum uncertainty (decision boundary)
191	0.75	High probability hypoxic
255	1.00	Certain hypoxic

Method	Equation	Range	Interpretation
Raw Probability	$P(Y=1) = \frac{1}{N}\sum_i T_i$	$[0, 1]$	Fraction of trees voting hypoxic
Calibrated Probability	$P_{\text{cal}} = \sigma(A \cdot f + B)$	$[0, 1]$	Adjusted for empirical frequencies
Shannon Entropy	$H = -p\log_2 p - (1-p)\log_2(1-p)$	$[0, 1]$ bits	Information content (0=certain, 1=uncertain)
Uint8 Encoding	$\text{pixel} = 255 \times p$	$[0, 255]$	Storage-efficient raster format

Field	Description
`file`	Date stamp from filename (e.g., `2020-06-01`)
`class_0_area_sqkm`	Non-hypoxic area (km²)
`class_1_area_sqkm`	Hypoxia zone extent (km²)

	file	class_0_area_sqkm	class_1_area_sqkm
0	2002-07-01	1.461363e+06	109060.962013
1	2002-08-01	1.487167e+06	83257.198658
2	2002-09-01	1.524132e+06	46291.944892
3	2002-10-01	1.511844e+06	58579.763991
4	2002-11-01	1.525554e+06	44870.078064
...	...	...	...
323	2009-03-01	1.536844e+06	33579.700187
324	2009-04-01	1.538245e+06	32179.177781
325	2009-05-01	1.537058e+06	33366.255975
326	2009-06-01	1.506612e+06	63812.430942
327	2009-07-01	1.441411e+06	129013.070202

Metric	Target	Current	Notes
Model F1 Score	>0.90	0.92	✅ Exceeds target
API Latency (p95)	<500ms	TBD	Measure after deployment
Dashboard Load Time	<3s	TBD	Test with 10+ rasters
Prediction Coverage	>90% pixels	TBD	Check cloud/nodata gaps
User Engagement	50 monthly users	0	Post-launch metric

Environment Set Up¶

Code and Data Download¶

Prerequisites¶

Create Virtual Environment¶

Step 1: Create Virtual Environment¶

Step 2: Activate Virtual Environment¶

Step 3: Verify Python Version¶

Step 4: Upgrade pip¶

Step 5: Install Required Packages¶

Step 6: Optional Geospatial Packages (uncomment in code if needed)¶

Step 7: Verify Installation¶

Table of Contents¶

Hypoxia Prediction in the Gulf of Mexico Using Machine Learning¶

Overview¶

Objectives¶

Data Sources¶

In-Situ Measurements¶

Satellite Imagery¶

Coordinate Reference Systems¶

Methodology¶

1. Data Loading & Exploration¶

2. Data Integration¶

3. Model Development¶

4. Model Interpretation & Uncertainty Quantification¶

5. Spatial Prediction on Rasters¶

6. Spatial Quantification & Temporal Analysis¶

Key Features¶

Requirements¶

Python Environment¶

Installation¶

Project Structure¶

Workflow Summary¶

Part 1: Data Loading & Exploration¶

Part 2: Data Integration¶

Part 3: Model Building¶

Part 4: Interpretation & Diagnostics¶

Part 5: Spatial Prediction¶

Part 6: Spatial Quantification¶

Part 7: Spatial-Temporal Trend Analysis¶

Hypoxia Prediction: Whole Pipeline Overview¶

Part 8: Next steps¶

Outputs¶

Model Artifacts¶

Prediction Rasters¶

Quantitative Results¶

Visualizations¶

Usage¶

1. Run the Notebook¶

Key Results¶

Model Performance (Test Set)¶

Feature Importance (Top 3)¶

Temporal Patterns¶

Validation & Limitations¶

Strengths¶

Limitations¶

Future Improvements¶

Citation¶

Contact¶

License¶

Part 1: Data Loading and Exploration¶

Exploratory Data Analysis¶

Do for filtering, adding hypoxia column, field data visualization¶

Checking decadal changce¶

Part 1. 2 : Extracting Satellite Data from Cloud and utilizing Cloud Computing¶

What We're Doing:¶

What do these variables mean?¶

Why This Matters for Hypoxia Prediction:¶

Why MODIS-Aqua L3SMI (NASA/OCEANDATA/MODIS-Aqua/L3SMI)?¶

Why Chlorophyll‑a, POC, and SST are sensitive to hypoxia¶

Generating Training points¶

Part 2: Data Integration¶

Part 3: Model Building & Training¶

Overview: Classical Machine Learning Approach (Teaching Purpose)¶

Related Published Work & Extensions¶

Current Implementation: Classical Baseline¶

Future work will explore adding one or more advanced models from the table above, depending on data availability and computational resources.¶

Why balance the classes?¶

Balancing approach used here¶

Train/test split rationale¶

Model registry & evaluation metrics¶

Code and Data Download ¶

Why MODIS-Aqua L3SMI (`NASA/OCEANDATA/MODIS-Aqua/L3SMI`)?¶

Partial Dependence Plot: Chlorophyll (`chlor_a`)¶

1. Probability Estimation via `predict_proba()`¶